TCGA BRCA
TCGA-BRCA data preparation
Download TCGA-BRAC data
# query <- GDCquery(project = "TCGA-BRCA",
# data.category = "Transcriptome Profiling",
# data.type = "Gene Expression Quantification")
#
query <- GDCquery(project = "TCGA-BRCA",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
sample.type = c("Primary Tumor","Solid Tissue Normal"))# unstranded: Strand 비특이적 발현 데이터
# stranded_first: Strand 특이적 (첫 번째 strand) 발현 데이터
# stranded_second: Strand 특이적 (두 번째 strand) 발현 데이터
# tpm_unstrand: TPM (Transcripts Per Million) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_unstrand: FPKM (Fragments Per Kilobase of transcript per Million mapped reads) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_uq_unstrand: FPKM-UQ (Upper Quartile normalized FPKM) 방식으로 정규화된 strand 비특이적 발현 데이터Select the raw reads
# ENSG to Symbols
library(org.Hs.eg.db)
# 유전자 ID에서 버전 제거
ensembl_ids <- gsub("\\..*", "", rownames(counts))
# 중복된 ID 처리
unique_ids <- unique(ensembl_ids)
# 중복 제거된 ID에 대해 유전자 심볼 매핑
gene_symbols <- mapIds(org.Hs.eg.db,
keys = unique_ids,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first")
# 원래 데이터에 매핑된 유전자 심볼 할당
# 여기서 중복 제거된 목록을 이용해 각각의 원본 ID에 매핑된 심볼 할당
symbol_names <- gene_symbols[ensembl_ids]
rownames(counts) <- symbol_namesTCGA BRCA 프로젝트의 메타데이터
제공된 열 이름 목록은 TCGA BRCA 프로젝트의 메타데이터에서 사용된 모든 열을 나타냅니다. 각 열의 용도는 다양하며, 환자와 샘플에 관련된 상세한 정보를 포함합니다. 이 정보는 연구자가 종양의 특성, 환자의 임상 경로, 생물학적 데이터와의 관계를 분석하는 데 사용됩니다. 아래는 각 열 이름의 의미와 사용 용도를 설명합니다:
- barcode: TCGA 프로젝트 내 각 샘플에 대한 고유 식별 코드입니다.
- patient: 각 환자에게 할당된 고유 식별자입니다.
- sample: 샘플의 고유 식별자입니다.
- shortLetterCode: 샘플의 간단한 코드를 나타냅니다.
- definition: 샘플 유형의 정의를 설명합니다.
- sample_submitter_id: 제출된 샘플의 ID입니다.
- sample_type_id: 샘플 유형의 ID입니다.
- tumor_descriptor: 종양의 설명을 제공합니다.
- sample_id: 샘플 ID입니다.
- sample_type: 샘플의 유형을 나타냅니다 (예: Primary Tumor).
- composition: 샘플 구성을 설명합니다.
- days_to_collection: 샘플 수집까지 걸린 일수입니다.
- state: 샘플의 상태를 나타냅니다.
- initial_weight: 샘플의 초기 무게입니다.
- preservation_method: 샘플 보존 방법입니다.
- pathology_report_uuid: 병리 보고서의 고유 식별자입니다.
- submitter_id: 제출자 ID입니다.
- oct_embedded: OCT(광학 단층 촬영)에 포함되었는지 여부입니다.
- specimen_type: 표본 유형입니다.
- is_ffpe: FFPE(형식고정 파라핀 포매) 처리 여부입니다.
- tissue_type: 조직 유형입니다.
- synchronous_malignancy: 동시 악성 종양의 유무입니다.
- ajcc_pathologic_stage: AJCC 병리학적 단계입니다.
- days_to_diagnosis: 진단까지 걸린 일수입니다.
- treatments: 처리 내용입니다.
- last_known_disease_status: 마지막으로 알려진 질병 상태입니다.
- tissue_or_organ_of_origin: 기원 조직 또는 기관입니다.
- days_to_last_follow_up: 마지막 추적 관찰까지 걸린 일수입니다.
- age_at_diagnosis: 진단 당시 나이입니다.
- primary_diagnosis: 주 진단입니다.
- prior_malignancy: 이전 악성 종양의 유무입니다.
- year_of_diagnosis: 진단 연도입니다.
- prior_treatment: 이전 치료 내용입니다.
- ajcc_staging_system_edition: AJCC 병기 시스템 에디션입니다.
- ajcc_pathologic_t: AJCC 병리학적 T 등급입니다.
- morphology: 형태학입니다.
- ajcc_pathologic_n: AJCC 병리학적 N 등급입니다.
- ajcc_pathologic_m: AJCC 병리학적 M 등급입니다.
- classification_of_tumor: 종양 분류입니다.
- diagnosis_id: 진단 ID입니다.
- icd_10_code: 국제 질병 분류 코드(ICD-10)입니다.
- site_of_resection_or_biopsy: 절제 또는 생검 부위입니다.
- tumor_grade: 종양 등급입니다.
- progression_or_recurrence: 진행 또는 재발 여부입니다.
- alcohol_history: 음주력입니다.
- exposure_id: 노출 ID입니다.
- race: 인종입니다.
- gender: 성별입니다.
- ethnicity: 민족입니다.
- vital_status: 생존 상태입니다.
- age_at_index: 인덱스 시 나이입니다.
- days_to_birth: 출생까지 걸린 일수입니다.
- year_of_birth: 출생 연도입니다.
- demographic_id: 인구통계학적 ID입니다.
- days_to_death: 사망까지 걸린 일수입니다.
- year_of_death: 사망 연도입니다.
- bcr_patient_barcode: BCR 환자 바코드입니다.
- primary_site: 주요 발생 부위입니다.
- project_id: 프로젝트 ID입니다.
- disease_type: 질병 유형입니다.
- name: 이름입니다.
- releasable: 공개 가능 여부입니다.
- released: 공개 여부입니다.
- days_to_sample_procurement: 샘플 조달까지 걸린 일수입니다.
- paper_patient: 논문에 사용된 환자 데이터입니다.
- paper_Tumor.Type: 논문에 기술된 종양 유형입니다.
- paper_Included_in_previous_marker_papers: 이전 마커 논문에 포함되었는지 여부입니다.
- paper_vital_status: 논문에 기술된 생존 상태입니다.
- paper_days_to_birth: 논문에 기술된 출생까지의 일수입니다.
- paper_days_to_death: 논문에 기술된 사망까지의 일수입니다.
- paper_days_to_last_followup: 논문에 기술된 마지막 추적 관찰까지의 일수입니다.
- paper_age_at_initial_pathologic_diagnosis: 논문에 기술된 최초 병리학적 진단 시 나이입니다.
- paper_pathologic_stage: 논문에 기술된 병리학적 단계입니다.
- paper_Tumor_Grade: 논문에 기술된 종양 등급입니다.
- paper_BRCA_Pathology: 논문에 기술된 BRCA 병리학입니다.
- paper_BRCA_Subtype_PAM50: 논문에 기술된 BRCA 하위 유형 PAM50입니다.
- paper_MSI_status: 논문에 기술된 MSI 상태입니다.
- paper_HPV_Status: 논문에 기술된 인유두종바이러스(HPV) 상태입니다.
- paper_tobacco_smoking_history: 논문에 기술된 흡연력입니다.
- paper_CNV Clusters: 논문에 기술된 CNV 클러스터입니다.
- paper_Mutation Clusters: 논문에 기술된 돌연변이 클러스터입니다.
- paper_DNA.Methylation Clusters: 논문에 기술된 DNA 메틸화 클러스터입니다.
- paper_mRNA Clusters: 논문에 기술된 mRNA 클러스터입니다.
- paper_miRNA Clusters: 논문에 기술된 miRNA 클러스터입니다.
- paper_lncRNA Clusters: 논문에 기술된 lncRNA 클러스터입니다.
- paper_Protein Clusters: 논문에 기술된 단백질 클러스터입니다.
- paper_PARADIGM Clusters: 논문에 기술된 PARADIGM 클러스터입니다.
- paper_Pan-Gyn Clusters: 논문에 기술된 Pan-Gyn 클러스터입니다.
TCGA-BRCA Analysis
Create DESeq object
PCA plot
ggplot(pcaData, aes(x = PC1, y = PC2, fill = sample_type)) +
geom_point(size = 2, alpha = 0.8, shape = 21, color = "black", stroke = 0.2) +
# ggrepel::geom_text_repel(aes(label=name),
# color="grey6", size=3, hjust= -0.3, vjust=-0.3) +
labs(x = paste("PC1: ", round(100 * PCA_var[1]), "% variance"),
y = paste("PC2: ", round(100 * PCA_var[2]), "% variance")) +
theme_bw() +
theme(legend.title = element_blank()) +
ggtitle("PCA") +
labs(caption = " ")Differentially Expressed Genes (DEG) analysis
DEG strategy
Differentially Expressed Genes (DEG) between Tumor tissue and Normal tissue
Tumor tissue/Normal tissue
n = 1224
Tumor n = 1111
Normal n = 113
Create DEG object
# Generate info table
info <- data.frame(matrix(nrow = ncol(count.mtx), ncol = 2))
colnames(info) <- c('sample', 'cond')
info$sample <- colnames(count.mtx)
info$cond <- dds$sample_type
info$cond <- factor(info$cond,
levels = c("Solid Tissue Normal","Primary Tumor")) # CTL going first
# levels(info$cond)
# DESeq
dds <- DESeqDataSetFromMatrix(count.mtx, info, ~ cond)
dds <- DESeq(dds)
# dds %>% saveRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)dds = readRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)# Add DEG information
fc = 1.5
pval = 0.05
# res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & pvalue < pval, 'UP',
# ifelse(log2FoldChange <= -log2(fc) & pvalue < pval, 'DN','no_sig')))
res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & padj < pval, 'UP',
ifelse(log2FoldChange <= -log2(fc) & padj < pval, 'DN','no_sig')))
res = na.omit(res)Volcanoplot
res$DE = factor(res$DE, levels = c('UP','DN','no_sig'))
res %>%
ggplot(aes(log2FoldChange, -log10(padj), color=DE)) +
geom_point(size=1, alpha=0.5) +
scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
theme_classic() +
geom_vline(xintercept = c(-log2(fc),log2(fc)), color='grey') +
geom_hline(yintercept = -log10(0.05),color='grey') +
guides(colour = guide_legend(override.aes = list(size=5))) +
ggtitle(paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )) +
ggeasy::easy_center_title() ## to center titleVolcanoplot with Number of DEGs
t= paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )
up = nrow(res[res$DE == "UP", ])
dn = nrow(res[res$DE == "DN", ])
res %>% ggplot(aes(log2FoldChange, -log10(padj), color=DE)) +
geom_point(size=0.5, shape=19, alpha=0.7) +
geom_vline(xintercept = c(-log2(fc), log2(fc)), size=0.1, color="grey") +
geom_hline(yintercept = -log10(0.05), size=0.1, color="grey") +
scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
theme_bw() +
annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up),
hjust = 1.1, vjust = 2, size = 5, color = "red") +
annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn),
hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
theme_bw() + ggtitle(t)Number of DEGs
res %>% filter(DE != "no_sig") %>%
ggplot(aes(DE, fill=DE)) + geom_bar(color="black", size=0.2) +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size= 4, color=c("salmon","royalblue")) +
scale_fill_manual(values = c("salmon", "royalblue"), guide=F) +
theme_bw()GSEA function
library(clusterProfiler)
hallmark <- msigdbr::msigdbr(species = "Homo sapiens", category = "H") %>%
dplyr::select(gs_name, gene_symbol)
perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
ranking <- function(res) {
df <- res$log2FoldChange
names(df) <- rownames(res)
df <- sort(df, decreasing = TRUE)
return(df)
}
ranked.res <- ranking(res)
x <- clusterProfiler::GSEA(geneList = ranked.res,
TERM2GENE = ref,
pvalueCutoff = pvalueCutoff,
pAdjustMethod = "BH",
verbose = TRUE,
seed = TRUE)
result <- x@result %>% arrange(desc(NES))
result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
return(result)
}
# Application
gsea.res = perform_GSEA(res = res, ref = hallmark) filtered_gsea = gsea.res %>% mutate(sig= ifelse(pvalue <= 0.05,"p value <= 0.05", "p value > 0.05"))
# Modified GSEA NES plot
gsea_nes_plot =function(gsea.res, title, fontsize.x = 5, fontsize.y = 6){
gsea.res %>% ggplot(aes(reorder(ID, NES), NES)) +
geom_col(aes(fill=sig), color="grey1", size=0.2) +
coord_flip() +
labs(x="Pathway", y="Normalized Enrichment Score",
title= "GSEA") +
theme_classic() +
# scale_fill_gradient(low = '#FF0000', high = '#E5E7E9') +
scale_fill_manual(values = c("#FF0000","grey88")) +
theme(axis.text.x= element_text(size=fontsize.x, face = 'bold'),
axis.text.y= element_text(size=fontsize.y, face = 'bold'),
axis.title =element_text(size=10)) +ggtitle(title)
}